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The  research  grant  was  initiated  on  July  1,  2003.  Since  then,  I  have  developed  the  theory  for  the  band 
structure  of  nanoscale  layered  semiconductor  structures  of  ZnO/MgO  and  GaN/AIN  with  the  inclusion  of 
Mn  impurities.  These  materials  have  the 

crystallographic  structure  of  wurtzites.  _ ^ 

While  it  has  been  usual  to  grow  these  structures 
along  the  c-axis  (the  growth  axis)  of  the  wurtzite 
structure,  this  leads  to  the  presence  of  polarization 
fields  that  are  generated  along  the  growth  axis  due  to 
internal  strain  and  due  to  the  spontaneous 
polarization  that  is  naturally  present  in  these 
structures.  This  internal  field  in  each  layer  affects 
the  band  edge  profile  of  the  layered  heterostructure. 

Recently  it  has  been  recognized  that  growth  of  Figure  1:  A  schematic  of  the  band  bending  in  layered 
planar  structures  along  the  R-plane  would,  in  structures  grown  along  the  c-axis,  due  to  spontaneous 

practical  terms,  remove  effects  of  the  polarization  polarization  and  piezo-electric  fields.  The  distorted 

fields  that  now  act  in  the  in-plane  direction.  This  bandedges  lead  to  poor  wavefunction  overlap, 
leaves  the  band  edge  profile  in  the  layered  structure 

to  be  rectangular.  The  presence  of  the  internal  electric  fields  for  the  case  of  c-axis  growth  leads  to  poorer 
wavefunction  overlap  and  therefore  such  structures  are  not  optimized  for  optical  applications. 


With  the  above  in  mind,  I  developed  the  theoretical  machinery  needed  to  account  for  the  band  structure  of 
the  layered  heterostructures  grown  along  any  crystallograpJiically  allowed  direction.  This  is  done  by 
performing  a  rotation  on  the  8x8  k.P  Lagrangian 
and  the  strain  Lagrangian.  The  material 
parameters  needed  for  the  two  are  known  and  we 
use  these  to  generate  the  band  structure  of 
layered  systems. 

The  solution  of  the  8  coupled  Schrodinger 
equations  in  the  k.P  theory  was  obtained  using 
the  finite  element  method.  This  method  has  been 
applied  extensively  by  me  over  the  past  decade 
to  implement  the  correct  boundary  conditions  in 
the  layered  III-V  heterostructures.  Now  this  is 
being  extended  to  the  wurtzites.  For  further 
details  on  the  application  of  finite  element 
methods  to  quantum  mechanics,  please  see  my 
book  published  by  Oxford  University  Press. 1  The 
finite  element  method  provides  the  flexibility  in  extracting  desired  levels  of  accuracy  for  energy  levels 
and  wavefunctions  through  the  number  of  finite  elements  used  in  each  layer.  It  also  permits  the  use  of  any 
boundary  condition  for  the  structure:  exponentially  falling  wavefunctions  for  quantum  wells,  traveling 
solutions  for  scattering  states  and  tunneling  boundary  conditions  for  open  structures. 

The  boundary  conditions  at  interfaces  require  the  maintenance  of  derivative  operator  ordering.  This  too  is 
implemented  in  the  framework  of  finite  elements.  I  have  developed  a  Lagrangian  approach  using  the 
action  integral  to  generate  Schrodinger’s  multiband  equations.  On  employing  Lowdin  perturbation  theory 
in  the  Lagrangian  formalism,  the  operator  ordering  issues  are  treated  in  a  very  natural  manner. 
Furthermore,  the  definition  of  the  conserved  current  and  the  current  continuity  condition  at  interfaces  can 
be  derived  unambiguously  by  employing  the  gauge-variational  approach  used  in  particle  physics  and  field 
theory. 
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Figure  2:  The  Wurtzite  heterostructures  grown  on  the 
R-plane,  i.e.,  with  layers  stacked  along  the  a-axis, 
show  no  spontaneous  polarization.  Thus  quantum 
wells  grown  along  this  direction  will  not  have  their 
bandedge  profiles  distorted  by  any  polarization  fields. 


In  parallel  with  this  effort,  I  have  been  developing  sparse  matrix  computational  methods  for  the  solution 
of  Schrodinger’s  equations  for  scattering/tunneling  structures  as  well  as  for  the  eigenvalue  problems  that 
arise  for  confined  states  in  quantum  well  structures.  We  have  fully  revised  our  codes  in  order  to  make 
them  compatible  with  ANSI-standards.  We  now  have  a  working  implementation  of  the  k.P  band  structure 
codes  for  wurtzite  structures. 

Polarization  effects:  Wurtzite  crystals  have  spontaneous  polarization  along  the  c-axis  and  also  inversion 
asymmetry  effects  represented  by  the  A7  valence  band  parameter.  This  parameter  is  known  to  be  large  in 
these  bulk  materials.  The  effect  of  this  inversion  asymmetry  is  to  generate  spin  splitting  between  the 
energy  doublets  which  exhibit  the  Kramers  degeneracy  in  the  absence  of  this  asymmetry.  The  energy 
bands  with  and  without  the  A7  term  are  shown  in  Fig.  3a.  When  the  heterostructure  is  grown  along  the 
[1120]  direction,  the  corresponding  rotated  Lagrangian  has  terms  proportional  to  A7  that  depend  on  kz’. 
The  wavevector  kz’  is  transformed  into  a  spatial  derivative  for  the  coordinate  along  this  new  growth  direct 
and  this  leads  to  contributions  to  the  energy  at  the  zone  center  in  quantum  wells.  This  is  the  origin  of  the 
energy  difference  at  the  zone  center  seen  in  Fig.  3b.  The  energy  bands  evolve  differently  in  this  case  as  is 
seen  on  comparing  the  figures  3a  and  3b.  These  effects  have  to  be  taken  into  account  in  designing  optical 
and  near-UV  lasers.  Furthermore,  they  could  be  important  when  we  start  modeling  terahertz  lasers  in  the 
wurtzite  materials  such  as  ZnO/MgZnO  (or  even  in  GaN/AlGaN  quantum  well  structures)  when  we 


1  L.  R.  Ram-Mohan,  “Finite  Element  and  Boundary  Element  Applications  to  Quantum  Mechanics,”  Oxford  University  Press  2002. 


consider  spintronics  effects.  We  have  been  able  to  model  this  effect  for  simple  quantum  wells  and 
superlattices. 

A  formal  paper  was  presented  at  the  International  Conference  on  the  Physics  of  Semiconductors  in  July 
(Publication  #4  below)  in  collaboration  with  workers  at  AFRL/WPAFB.  A  complete  report  on  this  work 
has  been  published  in  Superlattice  and  Microstructures. (#9  below). 

Magneto-Optics  in  DMS  structures:  I  have  been  working  towards  a  model  for  the  magneto-optical 
effects  in  ZnO./MgO  systems  to  build  in  the  effects  of  including  Mn  as  spin  dopants  in  the  layered 
structures.  Furthermore,  the  effect  of  free  carriers  and  the  approach  to  ferromagnetism  in  this  material 
system  will  be  modeled.  This  requires  the  presence  of  acceptors  in  specific  layers.  The  calculations  of 
energy  bands,  in  the  presence  of  a  magnetic  field,  allow  the  analysis  of  magneto-optical  data.  This 
together  with  the  Schrodinger-Poisson  selfconsistency  calculations  in  wurtzite  structures  is  necessary  to 
estimate  the  onset  of  ferromagnetism  in  heterostructures.  I  have  now  been  able  to  do  the  Schrodinger- 
Poisson  selfconsistency  in  an  arbitrary  layered  heterostructure  with  arbitrary  doping  profile  within  the 
framework  of  the  finite  element  method.  (#3,  in  the  list  of  paper  below).  The  DMS  interaction  and  the 
selfconsistent  band  bending  have  to  be  solved  simultaneously  in  order  to  explicitly  solve  for  the 
spontaneous  magnetization  in  DMS  systems.  We  have  now  been  able  to  complete  such  a  calculation  for 
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Figure  3:  (a)The  inplane  energy  dispersions  in  a  100  A  quantum  well  of  GaN/AlGaN  grown  along  the  [0001]  direction 
with  the  dispersion  curves  corresponding  to  the  presence  or  absence  of  the  inversion  asymmetry  term  A  7  in  the  valence 
bands  of  the  individual  layers.  Note  that  at  the  zone  center  the  spin  states  are  degenerate,  (b)  The  energy  dispersion  in  the 
case  of  a  100  A  GaN/AlGaN  quantum  well  grown  along  the  [1 120]  direction.  In  this  case  the  inversion  asymmetry  term 
induces  an  energy  shift  at  the  zone  center  but  the  spin  splitting  increases  as  the  wavevector  kx  increases. 


GAMnAs  quantum  wells  after  working  on  the  problem  for  full  two  year.  (#14,  below).  This  program  of 
work  will  be  continued  over  the  coming  year  even  though  the  grant  has  expired,  because  a  study  of  the 
onset  of  ferromagnetism  in  layered  DMS  heterostructures  is  one  of  the  key  problems  in  the  area 
spintronics. 


The  theory  of  LO-phonons  in  ternary  compounds  in  general,  and  wurtzite  materials  specifically,  was 
completely  reviewed  by  us  over  the  past  year.  We  found  that  the  energy  “dispersion”  as  a  function  of  the 
stoichiometric  composition  required  a  re-examination  because  of  the  one-mode  behavior  in  some 
ternaries  and  two-mode  behavior  in  others.  Most  of  the  earlier  work  insisted  on  using  the  two-mode 
treatment  with  the  random-element-iso-displacement  model,  and  this  led  to  the  mixed  modes  not  going  to 
the  appropriate  terminal  values  at  x=0  and  x=l  for  the  ion-displacement  amplitudes.  We  now  have 
sharpened  the  criteria  as  to  when  one  expects  a  one-mode  or  a  two-mode  behavior  in  these  materials 
which  models  to  use  when  describing  these  materials.  We  have  also  reexamined  the  Frohlich  electron 
phonon  interaction  to  obtain  the  correct  form  of  the  coupling  constant  in  the  case  of  ternary  compounds. 
In  effect,  this  re-examination  will  impact  on  some  of  the  modeling  that  has  already  been  reported  in  the 


•literature.  Tjiis  work  is  being  readied  for  publication,  as 
a  full  review  article  that  will  report  not  only  on  the  LO 
phonons  in  the  wurtzite  ternary  compounds  but  also  the 
zinc-blende  materials  and  systems  with  Mn-doping.  This 
work  is  still  in  progress. 

Magneto-Optics  in  Wurtzite  Structures:  I  am  now  in  a 

position  to  develop  the  full  theory  of  the  energy  levels  of 
electrons  in  layered  wurtzite  structures  in  the  presence  of 
a  magnetic  field.  I  distinguish  two  geometries:  in  the 
Faraday  geometry  the  magnetic  field  is  perpendicular  to 
the  layers,  and  in  the  Voigt  geometry  the  magnetic  field 
is  in  the  in-plane  direction.  The  latter  effect  has  not  been 
modeled  in  detail  in  heterostructures.  The  reason  is  that 
the  Landau  orbits  of  electronic  states  cross  the  interfaces 
and  there  is  no  simple  theory  available  that  can  give 
analytic  representations  for  the  wavefunctions.  I  have 
used  the  finite  element  method  to  model  this  effect  in  the 
zincblende  structures  over  the  past  year,  and  I  expect  to 
extend  this  theory  to  the  wurtzite  materials.  I  expect  to 
model  the  magneto-optical  effects  in  ZnO/MgO  systems 
and  build  in  the  effects  of  including  Mn  as  spin  dopants 
in  the  layered  structures.  Furthermore,  the  effect  of  free 
carriers  and  the  approach  to  ferromagnetism  in  this 
material  system  will  be  modeled.  This  requires  the 
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Fig.  4:  The  energy  dispersion  in  a  trilayer 
superlattice.  The  Finite  Element  approach  to 
modeling  the  band  structure  allows  for  any 
boundary  condition  at  the  ends  of  the 
heterostructure,  including  periodic  BC,  and  also 
allows  for  implementing  current  continuity 
conditions  at  internal  interfaces. 
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presence  of  acceptors  in  specific  layers.  The  calculations  of  energy  bands  in  the  presence  of  a  magnetic 
field  will  allow  analysis  of  magneto-optical  data.  This  is  an  integral  part  of  being  able  to  investigate 
ferromagnetic  behavior  in  spin-doped  heterostructures. 


As  a  final  note,  the  review  article  appearing  in  Journal  of  Physics  C  (#18  in  the  list  of  publications  below) 
was  published  in  November  2006.  It  provides  an  overview  of  the  work  done  here  and  details  of  the 
quantitative  theory  developed.  The  Editor  of  JoPC  sent  me  a  note  saying  that  over  just  the  month  that  this 
paper  has  been  available  there  have  been  250  downloads  of  copies  of  it.  According  to  him,  this  is  in  the 
top  10%  of  all  the  JoP  papers  in  all  fields  for  the  last  quarter  of  2006.  This  paper  is  attached  in  Appendix 
A  for  the  convenience  of  the  reader. 


In  all  of  the  above  work  I  have  had  active  interactions  with  Dr.  John  D.  Albrecht  who  is  at  the 
AFRL/WPAFB. 


Guidance  of  Personnel: 

1.  I  have  guided  the  work  of  Alexi  Girgis,  who  started  as  a  talented  undergraduate  student  at  WPI,  and  is 
now  a  graduate  student  with  me  since  2005.  I  gave  several  special  topics  courses,  over  and  beyond  my 
usual  course  load,  to  bring  him  up  to  the  level  where  he  could  do  research  with  me.  Alexi  collaborated  in 
this  work  and  has  provided  excellent  contributions  to  this  project  through  the  theory  development  and  the 
modeling  aspects.  Alexi  has  5  publications  to  his  credit  so  far,  which  is  an  excellent  record. 

2.  I  hired  a  Postdoctoral  Research  Associate,  Dr.  Rita  Gupta,  in  May,  2004,  and  her  contributions  will 
come  on  line  over  this  coming  year.  She  investigated  the  effect  of  phonons  in  the  wurtzite  structures  and 
their  effect  in  spin-polarized  cascade  lasers.  This  is  a  very  critical  aspect  of  modeling  the  onset  of 
ferromagnetism  in  the  wurtzite  spintronics  materials.  She  moved  to  the  University  of  Essex  in  August 
2005. 


3.  I  collaborated  with  Professor  K.  H.  Yoo  and  supported  his  sabbatical  visit  to  WPI,  for  work  on  this 

project.  Several  publications  have  come  out  of  this  collaboration.  We  are  planning  on  continuing  our 

collaboration  to  work  on  the  magneto-optics  of  layered  DMS  wurtzite  structures. 

Awards  and  Recognitions  received  during  this  grant: 

1.  Professor  Ram-Mohan  was  elected  to  the  Fellowship  of  the  Institute  of  Physics,  UK,  in  2003. 

2.  Ram-Mohan  was  awarded  the  Engineering  Excellence  Award  of  the  Optical  Society  of  America  for 
2003  “for  the  development  of  advanced  algorithms  for  the  successful  optimized  design  methods  for 
mid-IR  quantum  well  lasers,  for  advancing  the  paradigm  of  wavefunction  engineering,  and  for 
mentoring  students  in  optoelectronics.'1''  October  2003  at  the  OSA  annual  meeting  in  Tuscon,  AZ. 

3.  Ram-Mohan  was  awarded  the  Department  of  the  Air  Force  Certificate  of  Achievement,  given  “in 
recognition  of  the  outstanding  technical  accomplishments,  involving  fundamental  theoretical  issues  in 
optoelectronics  and  the  modeling  and  simulation  of  physical  systems,  made  during  2001-2003”  at  the 
Air  Force  Research  Laboratory,  Wright-Patterson  Air  Force  Base,  Dayton,  Ohio;  August  2003. 

4.  Ram-Mohan  was  awarded  the  Clark  Way  Harrison  Distinguished  Visiting  Professorship  at 
Washington  University  in  St.  Louis,  Missouri,  during  January-May,  2005. 

5.  Ram-Mohan  received  the  Sigma  Xi  Senior  Faculty  Award  for  Research  Excellence  at  WPI.  April 
2005,  “ given  in  recognition  of  his  development  of  the  finite  element  method  as  an  invaluable 
calculational  tool  in  the  fields  of  quantum  mechanics,  solid  state  and  semiconductor  physics,  and  for 
the  international  pre-eminence  that  he  has  earned  for  this  work." 

6.  Sun  Microsystems  grant:  Academic  Excellence  Award  of  an  equipment  grant  from  Sun  Microsystems 
of  $28,515  for  a  12  CPU  cluster  of  three  computers  for  my  work  on  wavefunction  engineering  of 
nanoscale  devices. 

7.  Student  Alexi  Girgis:  received  the  Sigma  Xi  Undergraduate  Research  Award  for  his  work  on  this 
project.  Alexi  is  continuing  towards  his  PhD  at  WPI  working  on  this  project. 

Publications  during  the  course  of  this  grant: 

1.  “Response  of  an  extraordinary  magnetoresistance  read-head  to  a  magnetic  bit,”  J.  Moussa,  L.  R.  Ram- 
Mohan,  A.  C.  H.  Rowe,  and  S.  A.  Solin,  Journal  of  Applied  Physics  94,  1110-1114  (2003). 
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Abstract 

We  present  the  theoretical  basis  of  wavefunction  engineering,  We  show  that 
a  Lagrangian  formulation  for  the  valence  bands  of  bulk  semiconductors  in 
the  k  •  P  model  provides  a  direct  approach  to  determining  derivative  operator 
ordering  in  the  multiband  description  of  electronic  states  in  semiconductors  in 
the  envelope  function  approximation.  The  current  continuity  condition  is  also 
obtained  through  a  gauge  variation  on  the  Lagrangian.  This  naturally  leads 
into  a  finite  element  approach  for  the  discretization  of  Schrodinger’s  equation. 
Furthermore,  by  including  the  Poisson  Lagrangian  a  self-consistent  treatment 
of  the  Schrodinger-Poisson  band-bending  in  arbitrarily  doped  structures  can 
be  calculated.  The  theory  is  developed  for  both  the  zinc  blende  and  wurtzite 
structured  compound  semiconductors  and  their  heterostructures.  Calculations 
for  quantum  wells  and  superlattices  are  presented  to  illustrate  wavefunction 
engineering  of  these  structures  and  the  control  achieved  in  obtaining  desirable 
wavefunction  localization.  We  show  that  when  combined  with  optimization 
methods,  wavefunction  engineering  provides  a  powerful  new  methodology  for 
the  optimized  design  of  optoelectronic  devices. 
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1.  Introduction 

There  have  been  remarkable  advances  over  the  past  few  decades  in  the  development  of 
optoelectronic  devices.  In  particular,  quantum  well  lasers,  emitting  with  wavelengths  from  the 
near-UV,  blue,  red,  the  mid-IR  region  (2-8  /tm),  all  the  way  to  the  far-IR  (~100  /tin)  region 
of  the  spectrum,  have  been  successfully  demonstrated.  This  spectacular  evolution  since  the 
development  of  GaAs  quantum  well  lasers  in  the  early  1970s  has  been  achieved  by  means  of 
technological  improvements  in  the  growth  of  structures.  In  the  following,  we  are  concerned 
with  the  electronic  properties  of  the  III — V  and  II — VI  compound  semiconductors  and  their 
heterostructures  [I], 

The  geometry  of  the  layered  structure,  i.e.,  the  layer  thicknesses  in  the  multi-quantum 
well  structure,  plays  a  key  role  in  determining  the  effective  energy  bandgap  of  the 
heterostructure,  which  arises  as  a  manifestation  of  the  Heisenberg  uncertainty  principle  and 
carrier  localization.  This  insight  was  distilled  into  the  concept  of  ‘bandgap  engineering’. 
By  choosing  different  materials  with  suitable  energy  bandgaps  for  individual  layers,  bandgap 
engineering  has  provided  a  conceptual  framework  for  designing  hctcrostructurcs  with  specific 
optical  properties  [2]. 

Through  a  combination  of  computational  modelling  and  advances  in  the  growth  of 
materials,  we  can  progress  considerably  further  today  in  exploring  the  details  of  the  carrier 
wavefunctions,  how  they  are  shaped  by  choices  for  materials  in  the  individual  layers,  and 
the  layer  geometry.  We  are  essentially  tailoring  the  localization  and  distribution  of  the 
wavefunctions  in  the  heterostructure,  and  optimizing  optical  transition  matrix  elements,  which 
fundamentally  determine  the  response  and  operation  of  the  optoelectronic  device.  In  this  article, 
we  provide  a  theoretical  and  computational  overview  of  these  advances  within  the  framework 
of  what  we  term  as  wavefunction  engineering  [3], 

We  have  emphasized  earlier  [4]  that  the  finite  element  method  (FEM)  can  be 
employed  to  advantage  in  simulating  the  electronic  energy  band  structures  of  semiconductor 
heterostructures.  The  FEM  [3,  5]  may  be  considered  to  be  the  discretization  of  the  action 
integral,  which  is  fundamental  to  all  of  physics  through  its  use  in  the  principle  of  least  action. 
We  show  that  the  Schrodinger  Lagrangian  (rather  than  the  Hamiltonian),  in  conjunction  with 
the  action  integral,  provides  specific  advantages  in  that  we  can  of  course  obtain  the  equation  of 
motion,  the  Schrodinger  equation.  We  can  advance  further  in  this  framework  by  employing  a 
gauge-variational  method  to  obtain  the  conserved  current  and  the  current  continuity  boundary 
conditions  at  interfaces  [6]  in  an  unambiguous  manner,  and  include  perturbations  such  as 
the  built-in  strain  in  the  hctcrostruclurc  and  externally  applied  electric  and  magnetic  fields 
in  a  natural  manner.  Again,  the  Lagrangian  may  be  used  in  a  discretized  action  in  FEM 
calculations  to  obtain  wavefunctions  and  momentum  matrix  elements.  In  the  layered  structure, 
each  layer  contributes  additively  to  the  global  action  integral.  Finally,  the  effect  of  charge 
redistributions  in  doped  heterostructures  due  to  ionized  impurity  dopants  and  free  carriers, 
leading  to  energy  bandedge  bending,  can  be  accounted  for  using  the  FEM.  This  ubiquitous 
problem  of  the  Schrodinger-Poisson  self-consistency  occurs  in  essentially  all  active  devices, 
and  the  FEM  provides  a  convergent,  stable  solution  in  tens  of  iterations  [7]  for  arbitrary  layered 
heleroslruclures  with  arbitrary  doping  profiles. 

When  combined  with  computational  optimization  routines  the  FEM  allows  us  to  predict 
the  material  properties  and  optimized  geometry  of  structures  under  bias  that  can  have  the 
desired  energy  level  structure  necessary  for  laser  emission  at  specific  wavelengths.  In  this 
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sense,  wavefunction  engineering  has  matured  to  a  level  where  we  can  predict  the  optoelectronic 
properties  of  a  given  heterostructure.  and  furthermore  provide  reliable  guidelines  for  the 
growth  of  heterostructures  having  specific  desired  outcomes  in  terms  of  their  optoelectronic 
response  [8]. 

This  review  is  organized  as  follows.  In  section  2,  we  describe  the  k  •  P  model  Lagrangian 
used  to  calculate  the  band  structure  of  heterostructures.  In  section  3,  we  show  how  the 
wavefunction  and  current  continuity  boundary  conditions  are  derived  consistent  with  operator 
ordering.  We  give  examples  of  wavefunction  engineering  in  section  4.  We  compare  the  present 
methodology  with  other  band  structure  modelling  in  section  5,  and  provide  concluding  remarks 
in  section  6. 


2.  A  Lagrangian  formulation  of  the  valence  band  structure 


For  our  purposes  here,  in  the  envelope  function  approximation  (EFA),  the  general  form  of 
the  wavefunction  may  be  considered  to  be  a  linear  combination  of  a  finite  number  of  band 
wavefunctions,  indexed  by  n,  of  the  form 

t(r)  =  f,(z)e'k'xe'k>yun(r )  =  F„(r)  u„(r).  (1) 


A  clear  perspective  on  the  formulation  of  the  band  structure  arises  on  using  the  principle  of  least 
action  in  setting  up  the  k  •  P  model  [9].  A  number  of  issues  such  as  derivative  operator-ordering 
are  resolved  in  a  very  natural  manner,  leading  to  a  number  of  insights  into  the  approximations 
that  are  invoked  in  the  modelling.  The  action  integral  that  leads  to  the  Schrodinger  equation  for 
the  valence  bands  is  given  by 


Hr  l2  i  pc 

d‘V(F»i«»)  3 —a  +  (V(r)-E)  (F„(r)u„(r))=  /dt/dVL. 

'  L  4m  (i  J  J  J 


(2) 


The  directed  derivatives  3  and  3  act  on  the  functions  appearing  to  the  left  and  to  the  right, 
respectively.  All  the  bands  are  treated  within  degenerate  perturbation  theory.  In  the  following, 
the  integral  over  time  is  irrelevant  since  we  are  concerned  with  the  time-independent  band 
structure  problem. 

Within  the  spirit  of  the  EFA,  we  perform  ‘cell-averaging'  by  integrating  over  each  unit 
cell  in  the  crystal.  The  envelope  functions  Fn(r)  are  typically  considered  to  be  slowly  varying 
functions  whereas  the  cell-periodic,  and  more  oscillatory,  Bloch  functions  satisfy  Schrodinger’s 
equation  with  band-edge  energies.  We  use  the  results: 


/ 
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and  write  the  momentum  matrix  element  between  Bloch  states  as 
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The  cell-averaged  action  integral  can  then  be  written  in  terins  of  the  envelope  functions  alone, 
and  we  have 


(A) 
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The  integrand  in  large  parentheses  in  the  above  equation,  the  Lagrangian  density  (L),  may  be 
cast  in  a  matrix  form  by  separating  the  band  index  n  into  the  valence  v,  and  the  energetically 
higher  remote  bands  with  an  index  r.  Following  Lowdin’s  perturbation  theory  [10],  we 
eliminate  the  remote  band  wavefunctions  in  favour  of  the  valence  bands  of  interest  by 
employing  the  equations  of  motion  for  the  remote  bands.  If  we  had  been  developing  the 
bulk  k  •  P  theory,  the  following  terms  would  correspond  to  the  Kane  model  [  1 ,  9]  for  the 
band  structure,  with  terms  quadratic  in  the  derivatives  (0(A2)  terms)  within  the  valence  bands. 
Ignoring  spin  for  the  moment,  we  have 


. . .  K 


3- — 3  +  (£,,  -  E)SVV. 
2  ni0 


/V 


(£  -  Er) 


Pry' 3  F u' , 


[u,  v'  =  X,  Y,  Z), 


(6) 


where  the  p-like  valence  Bloch  states  in  the  III— V  and  II— VI  semiconductors  have  been  denoted 
by  their  symmetry-specified  form  in  terms  of  X,  Y  and  Z.  We  have  neglected  terms  of 
higher  order  in  perturbation  theory  and  used  Lowdin's  approximations.  The  advantage  of  the 
above  derivation  is  that  it  explicitly  retains  the  order  of  the  derivative  operators,  each  term 
in  the  Lagrangian  is  Hermitian,  and  there  is  no  ambiguity  about  reordering  of  operators  for 
symmetrization.  The  focus  is  on  the  bands  in  the  immediate  vicinity  of  the  bandgaps  of  the  III — 
V  or  II— VI  materials.  We  now  follow  Kane’s  discussion  for  the  evaluation  of  the  contributions 
of  the  remote  bands  to  the  valence  bands.  Extensions  to  include  the  conduction  bands  follow 
schemes  used  earlier  by  explicitly  subtracting  the  contribution  of  the  conduction  band  states 
from  the  Luttinger  y  parameters  [II,  12], 


2. 1.  The  zinc  blende  semiconductors 

The  possible  intermediate  states  contributing  to  the  T|5  valence  states  belong  to  the 
energetically  higher  conduction  bands  of  symmetry  T,  (s-like  states),  F|2  states  of  symmetry 
{2 z2  -  .v2  -  y2,  x/3(x2  -  y2)),  T i5  states  of  the  form  (yz,  zx,xy),  and  the  r25  states.  These 
states  are  compatible  with  the  atomic  s,  p,  d,  and  f  orbital  contributions  to  the  valence  electronic 
structure  of  semiconductors,  respectively.  It  has  been  estimated  that  the  contribution  of  the 
f  orbitals  to  the  valence  electronic  bands  in  semiconductors  is  negligible  and  hence  the  r25 
intermediate  states  will  be  ignored.  Following  the  notation  of  Foreman  [13],  we  define 

r, 

a  =  -  (1/3/ho)  Y  K*lAvMI2/(£,.  -  Er). 

r 

I'l5 

Tt  =  -  (1/3/ho)  Y  l<*l/»>r)l2/(Er  -  Er), 

r 

and 

r,2 

5  =  -  (1/6»ho)  Y  W/Ll«r)|7(E„  -  Er).  (7) 

r 

We  then  determine  the  contributions  to  the  individual  matrix  elements  of  The  matrix  of 
derivative  operators  that  act  on  the  envelope  functions  [Fvx,  F„v,  Fvz )  can  be  shown  to  be  of 
the  form 
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in  the  envelope  function  basis.  Here  A  =  I  —  6<t  —  125,  5=1—  67T,  C]  =  65  —  6tr, 
and  C2  =  6n .  While  in  the  bulk  semiconductor  (where  the  envelope  functions  are  of 
the  form  exp(ik  •  r))  the  ordering  of  the  derivative  operators  is  redundant  [9,  14],  and 
the  way  in  which  the  band  parameters  are  grouped  is  irrelevant,  it  becomes  crucial  in  the 
heterostructure.  The  interface  boundary  conditions  depend  on  the  particular  ordering  of  the 
differentiation  operators  relative  to  the  band  parameters  that  vary  across  the  interface.  The 
parameters  <r,  5  and  it  are  simply  related  to  the  usual  Luttinger  parameters  [14]  and  we  have 
Y\  =  —  1  +  2rr  +  47r  4-  45,  y2  =  a  —it  +25.  and  y2  =  a  +it  —  5.  The  three  parameters  based 
on  the  symmetry  type  can  be  determined  from  experiments,  as  are  the  Luttinger  parameters,  so 
that  the  implementation  of  envelope  function  continuity  and  the  probability  current  continuity 
conditions  at  interfaces  now  becomes  feasible  [15,  16],  The  inclusion  of  the  electron  spin 
and  the  effects  of  spin-orbit  interaction  are  straightforward  and  we  obtain  a  6  x  6  valence 
Lagrangian.  The  spin-orbit  interaction  splits  the  six-fold  degeneracy  at  the  zone  centre  into 
a  four-fold  degenerate  heavy-hole  (hh)  and  light-hole  (Hi)  bands  of  T8  symmetry  with  total 
angular  momentum  J  =  3/2,  and  a  doubly  degenerate  split-off  (so)  band  of  Ty  symmetry.  The 


valence  Lagrangian  in  the  |  J, 

m  basis  is  given  by 
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where,  specializing  to  the  case  of  the  [001]  growth  direction,  we  have 


P  =  Yiq  +  9ty,a.  -  (£„  -  E)i  Q  =  Yikj  ~  2 nzy2dz, 

R  =  -73y2  (7  -  kj)  +  i2V3 y3  kxky\  S±  =  -j=  [*±C, 3,  +  d-C2k± 

=  37!  P;(C2  -2C,)*±  +  A±(C|  -2C2)3j  : 
c±  —  [*±(Ci  +  C2) +  3-(C|  +  C2)*±]  . 

Here  [yi,  y2,  y2,  C 1 ,  C2[  are  in  units  of  (/r/2«o)and  k±  =  kx±iky.  The  quantities  S±,  C±,  £± 
and  S±,  C±,  S±  differ  by  the  exchange  of  the  location  of  the  directed  derivatives  and  k- 
components  in  the  expressions.  Thus,  the  Hermiticity  of  the  Lagrangian  can  be  verified 
with  more  care  than  is  usually  needed.  All  the  parameters  introduced  above  are  expressible 
in  terms  of  the  three  quantities  a ,  5,  and  it ,  or  equivalently,  in  terms  of  the  Luttinger  y 
parameters.  In  the  bulk  semiconductor,  we  have  5+  =  — S+,  E_  =  ,/3/2 S_,  and  C±  =  0. 
Thus,  the  operator  ordering  leads  to  the  breaking  up  of  the  usual  parameters  into  component 
contributions  that  appear  to  contribute  differently  in  layered  systems  as  compared  with  bulk 
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semiconductors.  Additional  perturbations,  e.g.  the  strain  energy,  are  additively  included  in  the 
Lagrangian  [15,  16]. 


2.2.  The  wurtzite  semiconductors 

We  consider  the  bands  in  the  immediate  vicinity  of  the  bandgap  of  wurtzite  materials  at  the 
zone  centre.  The  valence  states  belong  to  the  [r5v  :  (X,  T)  +  T|V  :  { Z}}  representations.  The 
deeper  valence  bands  belong  to  the  F^.  :  [X2  —  K2,  2 XY}  and  the  T )v  :  { Z]  representations. 
The  nearest  higher  conduction  bands  belong  to  the  T.v  states  of  (X,  Y )  symmetry  and  the  Fjc 
states  transform  as  the  [Z]  representation.  Here  we  are  using  the  notation  of  [  1 7],  for  example. 

The  second-order  valence  band  terms  of  (6)  are  evaluated  as  done  by  Kane  for  bulk 
zinc  blende  semiconductors.  We  define  parameters  for  the  contributions  of  the  higher  F ;c 
conduction  bands  as  follows: 


o\ 


=  -L 


/  <r'  '  ,|/>,Y.y|r]")|2/(£„ 

m  o  p 


=  —  y"l(rjz||pz|r|sl>|2/(£,  -  Er,), 
m  a  n 


and 


<rt  =  —  r  1  l^x.Kir'S})(r!5,i/Jz|rjZ|)/(^I.  -  £r,). 

mn  “ 


(II) 


The  r?c  intermediate  conduction  states  do  not  contribute  to  the  valence  bands.  The  contribution 
of  the  r${XZ,  YZ\  intermediate  states  can  be  parameterized  by  writing 


x  z 

Y7. 


&i  =  —  Xj<r'  MpzirJ  '>r-/(£,-£rs); 
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xz\ 
YZ  I 


h  =  —  T  l<rj2,|px.r|F  ')|-/( Ev  -  Er,), 

»!<>  ^ 


r5 


and 


*3  =  —  £<r{Z| \px.  nr. 

mo  *—• 


|  A'Z 
YZ 


>(r,'  ’l/tzir,  ')/(£„-£  r5). 


(12) 


r5 


The  F(,  states  of  symmetry  {X2  —  T2,  2XT]  require  an  additional  parameter  to  represent  their 
contribution.  We  define 


X 


-Li<r^Vir 

'»()  r 


X2  -  Y2 
2  XY 


>\~/(.Er  -  £t6). 


(13) 


The  valence  band  Lagrangian  for  the  wurtzite  crystal,  in  the  (X,  Y,  Z)  basis,  is  now  given 
by  [18.  191 

9  A  L  j  9  x  +  9  v  M 1 9 


(L>  I’L"'  = 


+  i)  =M2il ;  +  (£,  -  F.) 

t)  yN  lit.v  +  A-  tv  [  ;>  v , 

!) -  AA  a  v  -f  it ,  AY  it  - , 


9  yA/|  9  A  +  9  v  L  |  9  v 
+  9zM29.  +  (£,.-  £) 

9-iV29v  +  9  v]V'9:, 


<).VN29,  +  9-N;9.v 

9v/V29  ;  +  9;/V^9v 

9  v  A/t  9  v  +  9  v  A/.)  9  v 
+  9  -/.29  -  +  (£„  —  £) 


(14) 
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where  the  parameters  £/,  M and  N:  are  given  by 

£i  =  I  —  2a \  —  4X;  L2  =  1  -  2<t2;  Nt  =  -2 a,  +  4k: 

N[  =  — 4A;  N2  =  -2cr3;  /V2  =  — 253;  (15) 

Mi  =  1  -  4A:  M2  =  I  —  2(5 1 ;  M3  =  1  -  2<S2, 

in  units  of  h2/2m0.  The  order  of  the  derivative  operators  with  /V[  and  N'2  terms  is  reversed  with 
respect  to  N\  and  N2  terms.  In  the  bulk  semiconductor,  the  operator  ordering  is  irrelevant  and 
we  obtain  the  bulk  parameters  N\  =  N\  +  N J  and  N2  =  N2  +  /VL  Furthermore,  as  has  been 
shown  by  Chuang  and  Chang  [20],  the  six-fold  symmetry  under  rotation  about  the  e-axis  leads 
to  the  relation  (L\  —  M\)  =  N\. 

The  valence  band  edge  is  not  triply  degenerate  in  wurtzite  semiconductors,  and  the  crystal 
field  splitting  leads  to  the  shifts  in  the  band  edge  energy  £,,  by  {A'|£(C.t\)|X)  =  (T| Z.(c  f  ,  |  Y)  — 
Ai.and  (Z]L(c.f.)|Z)  =0.  The  spin-orbit  splitting,  generated  by 

h 

Ln  =  7T^V'/  x  p  -<t 
4nigC- 

~  T(.so).vtfv  T  T(So)yOv  +  L(SO)rrr:,  (16) 

is  parameterized  by  the  relations  [20] 

(X|Z,(S0„|K)  =  — i  A2;  <T|L(so),|Z>  =  (Z|L(so)v|X)  =  -iA3.  (17) 

The  spinor  matrix  elements  of  the  standard  Pauli  spin  matrices  a  are  readily  included  in 
the  above  derivation.  With  these  definitions  in  place,  we  consider  the  inclusion  of  spin  in  the 
Lagrangian.  We  make  the  same  choice  as  in  [20]  for  the  valence  band  basis  states  with  spin.  The 
valence  band  Lagrangian,  including  the  crystal-field  splitting  and  the  spin-orbit  interaction,  is 
given  by 


where,  specializing  to  the  case  of  the  [0001]  growth  direction,  and  switching  the  sign  of  the 
energy  to  be  positive  for  the  valence  bands,  we  have 


P  =  (A,  +  A2  +  £,  —  £)+  or;  Q  =  (A,  —  A2  +  £,  —  £)  +  a, 
R  =  [M , 9-  +  A2A-jj]  +  (£„  -  £);  K  =  As(kx  +  i k,)2, 

S±  =  ~[B-N2k±  -  k±N$-\-,  S±  =  ~[d-N'2k±  -  k±N2d;], 

v  2  v  2 


a  =  j^9.  (A 


i  +  A3)9;  +  (A2  +  A4) 


A  =  v/2A3. 


Here  the  Bir-Pikus  parameters  [  I S]  are  identified  by  the  relations 


A  i  =  L2\ 


A2  =  My, 


A3  =  M2  —  L2;  A4  =  (£i  +  M|)/2  -  My, 
A5  =  (TVj  +/V|)/2;  A(,  =  (/V2  +  N'2)/ \pl. 


(19) 


(20) 


The  parameter  A(,  does  not  appear  in  the  Lagrangian  and  instead  we  have  to  contend  with  the 
determination  of  two  separate  parameters  (N2  and  N2)  in  order  to  proceed.  This  requires  that 
one  of  the  parameters  rr3  or  53  (see  (II)  and  (12))  be  known  from  ab  initio  calculations.  We 
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have  tentatively  used  S3  =  5) .  Mireles  and  Ulloa  [21]  have  one  parameter  less  in  their  analysis. 
The  quantities  S±  and  S±  differ  by  the  exchange  of  the  location  of  the  directed  derivatives  and 
/c-components  in  the  expressions. 

The  inversion  asymmetry  terms  [22]  in  the  vvide-gap  wurtzite  semiconductors  have  been 
shown  to  provide  an  improved  description  of  the  valence  bands.  The  valence  band  parameter 
Aj  has  been  expressed  in  the  convenient  form  rf7  =  (— ih/m0\/2)(X \px\Z)  by  Dugdale  et  al 
[23],  A  detailed  accounting  of  the  various  contributions  from  band  mixing  and  the  admixture 
of  d-orbitals  to  the  inversion  asymmetry  parameter  has  been  given  by  Lew  Yan  Voon  et  al  [  1 7], 
In  the  { XYZ }  orbital  basis  we  have 

/  0  0  \\flA-]kx  \ 

Linv=  0  0  k/2  A7LV  .  (21) 

\ — i \f2A-]kx  — i\/2  Ajky  0  J 

This  3x3  valence  band  matrix  is  doubled  to  accommodate  the  spin  degree  of  freedom  and  the 
resulting  6x6  Lagrangian  is  added  to  the  valence  band  Lagrangian  in  ( 18). 


3.  Interface  boundary  conditions  for  current  continuity 


The  action  integral  is  varied  with  respect  to  the  envelope  functions  F*  and  the  principle  of  least 
action  is  invoked  to  derive  the  six  coupled  Schrodinger  differential  equations  in  the  six-band 
zinc  blende  k-P  model.  Similar  considerations  hold  for  the  wurtzite  materials.  In  the  following, 
we  simplify  the  discussion  of  the  boundary  conditions  by  considering  just  the  valence  bands 
without  including  spin,  with  the  Lagrangian  given  by  (8).  In  a  layered  heterostructure,  each 
layer  contributes  a  piece  of  the  action  integral.  The  material  properties  change  across  layer 
interfaces  and  we  attach  a  layer  index  A  to  each  of  the  parameters  A,  B.  Cj,  and  Ci  appearing 
in  (8).  With  the  layer  growth  along  the  z-axis,  the  terms  in  the  Lagrangian  are  separated  into 
those  with  different  numbers  of  directed  derivatives  in  z.  The  action  integral  takes  the  form 

A  =  J  df  +  AzFi  +  (5,Q<X)  +  Qai+3-)  +  R(X)  +  (E]x)Sij  -  £)1  j  Fj. 

(22) 


Here  the  unit  matrix  is  denoted  by  1,  and  the  other  matrix  coefficients  in  the  Lagrangian  are 
given  by 
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(23) 


Observe  that  while  the  individual  linear  derivative  term  Q  is  not  Hermitian  by  itself,  the 
combination  of  the  two  terms  involving  Q  and  Q+  together  is  Hermitian.  One  of  the  boundary 
conditions  at  interfaces  is  the  continuity  of  the  envelope  wavefunction.  The  second  boundary 
condition  at  the  interfaces  is  the  continuity  of  the  probability  current  across  the  interfaces, 
which  we  evaluate  using  a  gauge-variational  approach  that  is  common  in  particle  physics. 
Following  Gell-Mann  and  Levy  [6],  we  substitute  Fj(z)  — >  elAI:) /•)(?)  in  the  Lagrangian  (22) 
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and  develop  a  variation  of  the  Lagrangian  with  respect  to  the  gauge  function  A(z).  In  the 
presence  of  operator  ordering,  the  conserved  current  is  obtained  in  the  form 


1 


*(L) 
h  V«(9;A) 


=  i  [-F/k'/’TF,  +  (dzF*)K^] Fj  -  F*Q^Fj  +  F/Q^+F,-  . 


(24) 


The  continuity  conditions  at  an  interface  at  z  =  0  are  then  given  by  F|0-  =  F|o+,  and 
KF'  +  QFlo-  =  KF'  +  QF|«+. 


4.  Examples  of  wavefunction  engineering  of  quantum  heterostructures 

The  discretized  Schrodinger  equation  is  obtained  by  the  FEM,  which  may  be  thought  of  as 
the  discretization  of  the  action  integral  [5],  The  physical  region  is  split  into  small  regions, 
Unite  elements,  in  each  layer,  and  the  wavefunctions  arc  expressed  in  terms  of  interpolation 
polynomials  with  as-yet  unknown  coefficients  in  each  finite  clement.  The  spatial  dependence 
is  integrated  out.  and  the  element  matrices  are  overlaid  in  a  global  matrix,  for  the  structure  as  a 
whole,  in  order  to  ensure  that  the  continuity  conditions  discussed  in  section  3  arc  implemented. 
We  minimize  the  action  by  varying  the  unknown  coefficients.  The  resulting  discretized 
Schrodinger  matrix  equation  is  solved  by  sparse  matrix  methods.  This  approach  ties  in  the 
action  to  the  numerical  work  and  provides  a  systematic  approach  to  improving  the  variational 
solution.  The  contributions  of  specific  layers  can  be  emphasized  in  FEM  by  discretizing  them 
appropriately.  We  have  applied  this  method  to  obtain  the  energy  eigenvalues  and  wavefunctions 
in  heterostructures,  in-plane  energy  dispersions,  and  optical  matrix  elements.  The  method  is 
also  extended  to  obtain  a  self-consistent  solution  to  the  Schrodinger-Poisson  problem  of  band 
bending  [7],  a  problem  of  particular  interest  in  the  context  of  spintronics  applications  where 
the  carrier  concentrations  in  GaMnAs,  for  example,  can  be  as  high  as  ~102n  cm-3  [24],  In  the 
following,  wc  consider  examples  of  finite  element  modelling  of  hctcrostructurcs  of  common 
interest.  The  calculations  make  use  of  the  band  parameters  given  in  [25], 

In  figure  I ,  the  evolution  of  the  design  of  a  heterostructure  for  a  mid-IR  laser  is  shown.  The 
initial  typc-11  quantum  well  structure  in  figure  1(a)  has  poor  conduction  band  (T'„)  and  valence 
band  (4'/,)  wavefunction  overlap  for  optical  matrix  elements.  By  introducing  another  layer  in 
the  well  to  form  a  step  potential,  as  in  figure  1(b),  the  wavefunctions  are  better  confined  and 
show  increased  overlap.  With  yet  another  layer,  we  obtain  the  so-called  W-structure,  shown  in 
figure  I  (c).  The  W-structure  was  grown  to  specifications  based  on  wavefunction  engineering, 
and  the  very  first  structure  grown  exhibited  lasing  in  the  mid-IR  region  (4  gm)  of  the  spectrum. 
This  is  a  classic  example  of  wavefunction  engineering  [4,  26]  through  which  we  are  able  to 
manipulate  the  material  properties  of  the  structure  to  obtain  desired  wavefunctions  and  design 
a  working  laser  to  operate  at  expected  wavelengths.  The  W-Iaser  holds  the  present  record  for 
continuous  operation  up  to  275  K  for  mid-IR  lasers. 

In  figure  2(a),  we  show  the  band  bending  in  a  modulation-doped  single  quantum  well, 
comparing  the  bandedge  profile  before  and  after  the  self-consistent  calculation,  with  the  rate 
of  convergence  for  the  Schrodinger-Poisson  calculation  shown  in  figure  2(b)  [7],  Figure  2(c) 
shows  the  first  heavy  hole  wavefunction  before  (dotted  lines)  and  alter  (solid  lines)  the  hole 
redistribution.  If  only  the  original  band  offset  potential  is  considered,  the  wavefunction  is  of 
the  simple  cosine-like  form.  When  the  self-consistent  potential  is  included,  the  total  potential 
resembles  that  of  a  double  quantum  well,  and  the  heavy  hole  wavefunction  reflects  this  by 
having  two  peaks.  In  more  complex  situations,  the  doping  profile  can  be  controlled  to  shape 
the  wavefunctions  differently,  a  situation  that  may  be  termed  wavefunction  engineering  by 
modulation  doping. 
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Type- f I  InAsf  1 7A)-CalnSb(3&A)  SL 

lnAs(24A)-GnInSb(31  A)-GnAISb(43A)  MOW 


lrAs(2 1  A)-GalnSb(31  A)-lnAs(?1  A)-GaA'SbM3A)  MQW 


(a) 


(b) 


(c) 


Figure  1.  (a)  The  lypc-ll  heterostruclurc  of  InAs/GalnSb  leads  to  localization  of  wavcfunclions 
in  different  layers.  This  is  ameliorated  partially  by  introducing  a  step  potential  as  in  (b).  The  W- 
struciurc,  shown  in  (c),  provides  a  substantially  improved  overlap  (figure  adapted  from  [4.  26);  with 
permission  from  Elsevier). 
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Figure  2.  The  valence  band  bending  in  a  quantum  well  with  p-doped  barriers  is  shown.  The 
FEM  iterations  for  die  Schrodinger-Poisson  calculations  converge  to  the  self-consistent  solution  in 
a  stable  manner  to  the  band  profile  shown  in  (a),  in  very  few  iterations  as  shown  in  (b),  in  the  dual¬ 
loop  algorithm  of  [7 1.  The  wavcfunclions  in  the  double  well  created  by  the  self-consistent  band 
bending  are  shown  in  (c)  (after  [7p. 


The  ZnSe/Zni_j.Mn,vSe  multilayer  system  was  reported  [27]  to  form  a  spin  superlattice. 
For  x  ~  0.04.  the  conduction  bandedge  of  ZnSe  is  almost  the  same  as  that  of  the  diluted 
magnetic  semiconductor  (DMS)  Zni_vMnvSe.  In  an  applied  magnetic  field,  the  Mn  spins 
are  aligned  and  the  effective  conduction  bandedge  of  ZnMnSe  for  spin-up  electrons  is  higher 
than  that  of  ZnSe.  It  is  the  opposite  for  spin-down  electrons,  so  that  the  spin-up  electrons  are 
localized  in  the  ZnSe  layer  and  spin-down  electrons  in  the  ZnMnSe  layer.  In  figure  3(a),  we 
plot  the  effective  bandedge  and  the  wavefunction  for  the  lowest  spin-down  electron  Landau 
level,  and  in  figure  3(b),  for  the  lowest  spin-up  level,  at  B  =  5  T. 

When  the  Mn  alignment  is  not  considered,  by  neglecting  the  DMS  Lagrangian  for  the 
DMS  layer  (dotted  lines),  the  probability  that  the  electrons  reside  in  the  ZnSe  layer  is  51.4% 
for  both  spin-up  and  spin-down  electrons  due  to  a  small  conduction  band  offset.  When  the  Mn 
alignment  is  taken  into  account  by  including  the  DMS  Lagrangian  for  the  ZnMnSe  layer,  the 
probability  becomes  42. 1  %  for  spin-down  electrons  and  60.6%  for  spin-up  electrons.  Thus  the 
Mn  alignment  through  an  external  magnetic  field  alters  the  localization  of  spins  in  different 
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Figure  3.  The  conduction  bandedge  profile  without  (dotted  lines)  and  with  (solid  lines)  Mn  spin 
alignment  in  a  magnetic  field  of  5  T  for  (a)  the  spin-down  and  (b)  for  the  spin-up  electron  in  a 
ZnMn/ZnMnSe  superlattice  is  shown.  The  corresponding  wavefunctions  arc  shown  in  the  *spin- 
supcrlattice’  formed  in  the  structure. 
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Figure  4.  The  magnetization  due  to  Mn  spin  alignment  in  a  100  A  CiaMnAs  quantum  well  is  shown. 
In  (a),  the  alteration  of  magnetization  with  increasing  hole  concentrations  from  5  x  I017  cm-3- 
10|l)  cm-3  is  shown.  With  further  increases  in  hole  concentration,  (b)  and  (c),  the  magnetization 
oscillates  in  orientation  across  the  quantum  well.  This  is  due  to  the  occupied  energy  levels  and  the 
corresponding  wavefunctions  for  holes  in  the  self-consistent  potential  (after  [24]). 


layers  in  specially  designed  DMS  superlattices.  This  is  an  example  of  magnetically  tuned 
wavefunction  engineering. 

Figure  4  shows  the  spontaneous  magnetization  of  Mn  ions  in  a  100  A 
Guo.t15Mno.05As/Alo.35Gao.65As  quantum  well  as  a  function  of  hole  concentration,  which  can 
vary  independently  of  the  M11  concentration.  Since  the  exchange  interaction  between  Mn  ions 
and  itinerant  holes  increases  with  hole  concentration,  the  magnetization  increases  with  hole 
concentration,  and  saturates  with  a  carrier  density  at  about  3  x  1019  cm”3.  The  shape  of  the 
magnetization  is  a  reflection  of  the  difference  between  spin-up  and  spin-down  hole  densities. 
As  the  hole  concentration  increases,  higher  subbands  are  occupied,  and  the  shape  of  magne¬ 
tization  changes  as  the  relative  distribution  of  spin-up  and  spin-down  holes  changes.  These 
Schrodinger-Poisson-DMS  self-consistent  calculations  were  performed  in  the  valence  four- 
band  Luttinger  model  using  the  FEM,  and  illustrate  the  power  of  the  method  in  successfully 
accounting  for  such  large  carrier  concentrations  [24]. 

We  now  consider  an  example  of  wavefunction  engineering  in  a  1 00  A  wurtzite  GaN/AlGaN 
quantum  well.  When  the  structure  is  grown  along  the  c-axis  the  resulting  polarization  leads 
to  a  distortion  of  the  quantum  well  profile  due  to  the  large  internal  piezoelectric  polarization 
fields.  The  wavefunctions  are  then  displaced  in  position  and  the  wavefunction  overlap  between 
electrons  and  boles  is  substantially  reduced,  as  shown  in  figure  5(a).  When  the  growth  is 
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Growth  Direction  [0001]  (A)  Growth  Direction  [1120]  (A) 

(a)  (b) 

Figure  5.  The  remarkably  high  internal  piezoelectric  fields  in  wurtzitc  GaN/AIGaN  structures 
distort  the  band-edge  potential  in  quantum  wells,  as  shown  in  (a),  leading  to  a  poor  overlap  of 
carrier  wavefunctions.  When  grown  along  the  [1 120]  direction  the  piezoelectric  fields  are  along  the 
plane,  leading  to  belter  characteristics  for  quantum  well  lasers.  Strain,  piezoelectric  effects,  and 
inversion  asymmetry  effects  have  been  included  in  the  calculations  (after  f  10].  Figure  reproduced 
with  permission  from  Elsevier). 


along  the  [1 1201  direction,  the  spontaneous  and  piezoelectric  polarizations  are  directed  in  the 
plane  (still  along  the  c-axis.  which  is  in  the  plane).  The  quantum  well  potential  is  no  longer 
distorted  by  the  internal  electric  field  and  the  wavefunction  overlap  is  then  ‘engineered’  to  be 
large.  The  calculation  of  the  energy  levels  and  wavefunctions  proceeds  by  employing  a  rotated 
Lagrangian  with  the  new  c'-axis  directed  towards  the  corresponding  growth  direction.  The 
coordinate  transformation  x'  =  Rx  leads  to  the  new  rotated  Lagrangian  given  by 

L|V..v'..-l(k,>  =  RL|.v.v.Z)(R_lk')R_l.  (25) 

The  matrix  elements  of  the  transformed  Lagrangian  in  the  envelope  function  basis  are 
separated,  using  symbolic  algebra  software,  into  terms  involving  powers  of  k[  =  -ifl/i);' 
corresponding  to  the  new  growth  direction.  The  rotated  Lagrangian  is  then  used  in  the  finite 
element  calculations  [18],  A  discussion  of  the  rotation  of  the  k  •  P  Hamiltonian  was  given 
earlier  in  fl  I.  28,  29,  21],  This  feature  of  the  absence  of  bandedge  profile  distortion  for  the 
[1120]  growth  direction  in  GaN/AIGaN  quantum  wells  has  been  investigated  experimentally 
by  [30,  31], 

As  a  final  example,  we  consider  the  design  of  an  intersubband  quantum  cascade  laser 
(QCL)  structure  for  emitting  in  the  terahertz  (THz)  region  of  the  spectrum.  Lasers  in  this 
region  of  the  spectrum  have  special  problems  to  overcome.  The  electronic  transitions  for 
photon  emission  are  in  the  energy  range  5-20  meV,  so  that  the  system  has  to  be  cooled  to 
cryogenic  temperatures,  and  the  laser  cavity  has  to  be  of  the  order  of  several  hundred  microns. 
In  designing  a  quantum  cascade  laser,  the  goal  of  the  optimization  is  to  have  a  three-level 
system  in  which  the  carriers  undergo  a  transition  from  energy  E2  to  E2,  emitting  a  photon.  The 
electrons  are  then  depleted  from  E2  by  a  transition  from  E2  to  Et  on  emitting  interface  and 
layer-confined  longitudinal-optic  phonons.  The  carriers  then  tunnel  into  the  next  period  of  the 
cascade,  arriving  in  the  uppermost  level  3  of  the  next  stage.  Recently,  the  operation  of  THz 
lasers  has  been  successfully  demonstrated  [32-34],  Here  we  are  concerned  with  the  optimized 
design  of  a  single  period  of  the  cascade  laser. 

In  order  to  design  such  a  structure  wc  consider  a  single  period  of  a  five-quantum  well 
structure  with  the  requirements  that  under  an  externally  applied  bias  of  say  2.5  kV  cm-1 
the  energy  differences  are  16  meV  (~77  /no  or  ~3.9  THz)  for  the  photonic  transition,  and 
~38  meV  for  the  phonon  emission.  The  calculation  proceeds  by  coupling  the  one-band 
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Figure  6.  A  five-well  structure,  representing  one  period  of  a  QCL  with  energy  level  structure  for 
emitting  photons  with  energy  16  meV  (3.9  THz)  is  shown  before  layer  thickness  optimization  in  (a). 
In  (b)  the  structure  after  GA  optimization  is  shown  with  well  and  barrier  widths  adjusted  to  obtain 
the  energy  levels  under  bias  that  arc  mentioned  in  the  text. 


finite  element  calculations  of  the  energy  levels  of  the  nine-layer  period  to  a  Powell  or  genetic 
algorithm  (GA)  optimizer  [35].  In  figure  6(a)  we  show  the  initial  well  widths  and  energy  levels 
with  no  bias.  Initially,  the  well  and  barrier  widths  are  uniformly  the  same  for  each  well,  and  the 
wavefunctions  have  been  shifted  by  the  eigenvalues  for  purposes  of  this  display.  This  may  be 
compared  with  the  final  structure,  figure  6(b).  obtained  through  GA  optimization  for  the  well 
and  barrier  widths  in  order  to  obtain  the  desired  energy  level  structure  in  GaAs/Alo.jGaojAs. 
The  wavefunctions  have  been  shifted  vertically  to  have  their  ,v-axes  coincide  with  the  energy 
levels  for  display  purposes.  The  gain  in  the  QCL  structure  per  period  was  not  optimized  with 
respect  to  the  rate  equations,  and  this  work  is  as  yet  preliminary  [36,  37]  in  this  regard. 

We  note  that  the  use  of  LO-phonon  emission  for  depleting  the  carriers  from  to  E\  can 
also  be  optimized  through  the  use  of  interface  phonon  modes  and  requiring  that  (En  -  E\) 
lies  in  the  set  of  interface  mode  energies.  The  interface  modes  extend  over  the  entire  structure 
and  hence  provide  better  overlap  with  the  carriers.  In  this  sense,  we  are  employing  phonon 
wavefunction  engineering  in  the  operation  of  the  QCL  [37], 

5.  Comparison  with  other  methods 

We  have  focused  mainly  on  the  k  •  P  model  and  the  use  of  the  envelope  function  modelling  for 
the  electronic  band  structure  near  the  centre  of  the  Brillouin  zone  (BZ).  Here  we  comment  on 
the  oilier  alternatives  for  evaluating  the  electronic  states  in  layered  heterostructures  so  as  to  put 
their  merits  in  perspective. 

(a)  In  the  empirical  tight  binding  model  (ETBM),  originally  put  forward  by  Slater  and 
Koster  [38],  the  electronic  states  are  considered  to  be  linear  combinations  of  atomic  (s, 
p,  d,  . . .)  orbitals.  The  Hamiltonian's  matrix  elements  between  the  atomic  orbital  states 
are  not  evaluated  directly,  but  are  instead  introduced  as  free  parameters  to  be  determined 
by  fitting  the  band  gaps  and  band  curvatures  (effective  masses)  at  critical  points  in  the  BZ. 
Depending  on  the  number  of  orbitals  and  nearest  neighbours  used  to  represent  the  states, 
the  ETBM  requires  that  the  overlap  integrals  be  determined  in  terms  of  the  measured  direct 
and  indirect  band  gaps  and/or  effective  masses  in  the  bulk  material  [39].  For  example,  the 
spV  basis  with  the  second-nearest-neighbour  scheme  turns  out  to  have  27  parameters 
for  the  zinc  blende  lattice  structure,  and  the  energies  and  effective  masses  are  obtained 
from  the  diagonalization  of  the  Hamiltonian.  It  has  been  emphasized  by  Lew  Yan  Voon 
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and  Ram-Mohan  [40]  that  the  bulk  band  curvatures  must  be  reproduced  well  in  order  to 
extend  the  analysis  reliably  to  heterostructures  since  the  effective  masses  play  a  crucial 
role  in  determining  the  energy  levels  in  quantum  wells,  for  example.  Thus,  reproducing 
the  bulk  bandgaps  alone  is  not  sufficient.  The  bulk  energy  bands  are  nonlinear  functions 
of  the  ETBM  parameters,  which  can  be  lilted  by  trial  and  error  or  by  using  optimization 
methods,  such  as  genetic  algorithms  [41], 

The  lack  of  a  direct  relationship  between  the  input  parameters  and  the  experimentally 
determined  quantities  is  probably  the  single  greatest  disadvantage  of  the  tight-binding 
method  in  making  complicated  band  structure  calculations.  The  strategy  for  employing  the 
ETBM  for  heterostructures  is  analogous  to  the  FEM  in  that  we  overlay  atomic  coupling 
matrices  (the  sub-Hamiltonians)  in  order  to  produce  a  global  matrix  for  the  heterostructure. 
The  ETBM  has  been  applied  to  superlattice  band  structure  calculations  [42],  A 
superperiodicity  boundary  condition  is  imposed  on  the  global  matrix,  and  the  interface 
boundary  conditions  correspond  to  interpolated  overlap  integral  parameters.  It  was 
originally  thought  that  the  momentum  matrix  elements  for  optical  transitions  require  a 
refitting  of  the  parameters  for  these  matrix  elements  as  well;  however,  it  was  shown  by  Lew 
Yan  Voon  and  Ram-Mohan  that  one  can  obtain  the  optical  matrix  elements  by  invoking 
the  Feynman-Hellmann  theorem.  We  can  relate  the  optical  momentum  matrix  element 
to  the  derivative  of  the  ETBM  Hamiltonian  with  respect  to  the  wavevector  so  that  the 
same  set  of  parameters  that  define  the  Hamiltonian  can  be  reused  in  evaluating  the  optical 
matrix  elements  in  a  self-consistent  manner  [40],  The  ETBM  provides  contributions  from 
the  entire  BZ  to  the  heterostructure  energy  levels,  whereas  the  k  ■  P  method  includes 
contributions  only  near  the  centre  of  the  BZ.  The  computational  requirements  for  2D 
and  3D  modelling  within  the  ETBM  would  be  more  demanding  in  our  opinion  than  the 
envelope  function  approach  with  FEM.  Inclusion  of  electric  and  magnetic  field  effects  arc 
also  much  more  complex  in  the  ETBM  than  in  the  procedure  we  have  reviewed  in  this 
article. 

(b)  While  the  inclusion  of  additional  bands  and  overlaps  of  higher  orbitals  is  a 
possible  approach  to  improving  the  tight-binding  modelling  of  energy  bands  in  bulk 
semiconductors,  the  effective  bond  orbital  model  (EBOM)  [43]  uses  spin-doubled  s,  pv, 
pv.  Pi  orbitals  to  generate  an  8  x  8  Hamiltonian.  A  crucial  difference  between  the  EBOM 
and  related  tight-binding  formulations  is  that  the  s  and  p  orbitals  are  centred  on  the  face- 
centred  cubic  lattice  sites  of  the  zinc  blende  crystal  rather  than  on  both  of  the  two  real 
atoms  per  lattice  site.  The  resulting  somewhat  ad  hoc  formulation  offers  considerable 
computational  savings  in  comparison  with  the  ETBM.  However,  the  main  significance  of 
the  EBOM  approach  derives  from  the  fact  that  the  resulting  secular  matrix  has  a  small- 
k  expansion  that  exactly  reproduces  the  form  of  the  eight-band  k  •  P  Hamiltonian.  This 
allows  the  EBOM  input  parameters  to  be  readily  expressed  in  terms  of  the  experimentally 
measured  parameters,  such  as  the  band  gap,  the  split-off  gap,  and  the  zone-centre  mass  of 
each  band,  which  has  not  been  accomplished  using  the  more  involved  ETBM.  In  fact,  the 
EBOM  can  be  thought  of  as  an  extension  of  the  k  •  P  method  to  provide  an  approximate 
representation  of  the  energy  bands  over  the  full  BZ.  Since  short-period  superlattice  bands 
sample  wavevectors  throughout  the  BZ,  we  may  expect  the  EBOM  to  be  more  accurate 
than  the  k  •  P  model  for  thin-layer  structures.  However,  the  EBOM  is  considerably  less 
efficient  computationally  than  the  k  •  P  method,  especially  for  thicker  superlattices.  Each 
lattice  position  must  be  represented  in  the  supercell  technique,  i.e.,  no  envelope  function 
approximation  is  made. 

(c)  The  influence  of  core  electrons  in  keeping  the  valence  electrons  outside  of  the  core  may  be 
represented  by  an  effective  repulsive  potential  in  the  core  region.  When  this  is  added  to  the 
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attractive  ionic  potential,  the  net  ‘pseudopotentiaF  nearly  cancels  at  short  distances  [44], 
The  valence  states  are  orthogonal  to  the  core  states,  and  the  resulting  band  structure  theory 
corresponds  to  the  nearly  free-electron  model.  In  the  empirical,  pseudopotential  model ,  the 
crystal  potential  is  represented  by  a  linear  superposition  of  atomic  potentials,  which  are 
modified  to  obtain  good  fits  to  the  experimental  direct  and  indirect  band  gaps  and  effective 
masses.  Further  details  are  presented  by  Cohen  and  Chelikowsky  [45]  and  in  the  reviews 
by  Heine  and  Cohen  [46],  Ab  initio  approaches  employ  calculated  band  parameters  (from 
the  density-functional  theory)  in  lieu  of  experimental  data.  Combinations  of  ab  initio  and 
empirical  methods  have  been  developed  to  a  high  level  of  sophistication  [47],  Extension 
of  the  pseudopotential  method  to  heterostructures  entails  the  construction  of  a  supercell  to 
assure  the  proper  periodic  boundary  conditions.  With  atomic  potentials  as  the  essential 
input,  the  electronic  properties  of  the  heterostructure  can  be  determined,  although  the 
required  computational  effort  far  exceeds  the  demands  of  the  k  •  P  method.  The  relative 
merits  of  the  k  •  P  and  pseudopotential  approaches  have  been  assessed  [48],  Again,  the 
inclusion  of  states  from  the  entire  BZ  in  the  construction  of  heterostructure  quantum  states 
in  these  schemes  entails  heavy  computational  resources. 


6.  Concluding  remarks 

We  have  shown  that  the  principle  of  least  action  can  be  employed  to  advantage  in  the  modelling 
of  semiconductor  heterostructures  by  the  FEM.  This  provides  a  natural  way  of  including 
complex  boundary  conditions  into  the  simulation.  The  FEM  is  not  a  low-accuracy  method, 
but  rather  it  can  be  used  to  solve  quantum  mechanical  problems  with  double  precision  accuracy 
at  the  nanoscale.  This  modelling  approach,  implemented  by  us  over  the  past  two  decades,  now 
allows  us  to  explore  optoelectronic  properties  of  complex  heterostructures  through  numerical 
simulations  [4,  49],  The  wavefunction  overlap,  optical  matrix  elements,  strain-induced  band 
splittings  leading  to  wavefunction  mixing,  energy  level  and  wavefunction  manipulation  through 
modulation  doping,  use  of  magnetic  impurity  layers  to  alter  distribution  of  the  wavefunctions, 
and  the  application  of  external  perturbations  such  as  electric  or  magnetic  fields  are  all 
increasingly  of  interest  in  controlling  the  actual  shape  of  the  wavefunctions.  These  represent 
aspects  of  wavefunction  engineering.  The  properties  of  the  structure  can  now  be  simulated  and 
optimized  on  the  computer  before  a  heterostructure  is  actually  grown  [  1 2,  26],  The  fundamental 
shift  in  paradigm  represented  by  wavefunction  engineering  is  leading  to  the  exploration  of  new 
electronic  mechanisms  and  concepts,  the  exploration  of  basic  physics,  and  a  rapid  turn-around 
in  the  design-growth-characterization  cycle  for  developing  new  optoelectronic  devices  with 
2D  (quantum  wire)  and  3D  (quantum  dot)  confinement  of  carriers.  We  are  just  beginning  to 
explore  the  freedom  in  what  may  be  called  heterostructure  architecture,  in  the  design  of  new 
optoelectronic  devices  in  three  dimensions  through  this  combination  of  fundamental  theory  and 
FEM  computations.  Such  an  approach  is  indispensable  and  perhaps  inevitable  for  modelling 
nanoscale  systems. 
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